%%Code for Element Generation

function elout=elcalc(elin,xn,yn,mat)
err=1e-12;
elout=elin;
elout(:,4)=1/2*(xn(elout(:,2)).*yn(elout(:,3))+xn(elout(:,1)).*yn(elout(:,2))+...
    xn(elout(:,3)).*yn(elout(:,2))+...
    xn(elout(:,3)).*yn(elout(:,1))-xn(elout(:,2)).*yn(elout(:,1))-...
    xn(elout(:,3)).*yn(elout(:,2))-xn(elout(:,1)).*yn(elout(:,3)));
izero=find(abs(elout(:,4))<err);
elout(izero,:)=[];
[ineg,jneg]=find(elout(:,4)<0);
elout(ineg,1:2)=[elout(ineg,2) elout(ineg,1)];
elout(ineg,4)=-elout(ineg,4);
elout(:,5)=sqrt((xn(elout(:,1))-xn(elout(:,2))).^2+...
    (yn(elout(:,1))-yn(elout(:,2))).^2);
elout(:,6)=sqrt((xn(elout(:,2))-xn(elout(:,3))).^2+...
    (yn(elout(:,2))-yn(elout(:,3))).^2);
elout(:,7)=sqrt((xn(elout(:,3))-xn(elout(:,1))).^2+...
    (yn(elout(:,3))-yn(elout(:,1))).^2);
elout(:,8)=(max([elout(:,5) elout(:,6) elout(:,7)]')./...
    min([elout(:,5) elout(:,6) elout(:,7)]'))';
den=(yn(elout(:,1)).*(xn(elout(:,3))-xn(elout(:,2)))+...
    yn(elout(:,2)).*(xn(elout(:,1))-xn(elout(:,3)))+...
    yn(elout(:,3)).*(xn(elout(:,2))-xn(elout(:,1))));
ij=find(den==0);
if(~isempty(ij))
    figure(1)
    hold on
    plot(xn(elout(:,1)),yn(elout(:,1)),'x')
    plot(xn(elout(:,1)),yn(elout(:,1)),'o')
    plot(xn(elout(:,2)),yn(elout(:,2)),'x')
    plot(xn(elout(:,2)),yn(elout(:,2)),'o')
    plot(xn(elout(:,3)),yn(elout(:,3)),'x')
    plot(xn(elout(:,3)),yn(elout(:,3)),'x')
    plot(xn(elout(:,3)),yn(elout(:,3)),'o')
    disp('Halted')
    pause
end
xcnum=(xn(elout(:,1)).^2+yn(elout(:,1)).^2).*(yn(elout(:,2))-...
    yn(elout(:,3)))+...
    (xn(elout(:,2)).^2+yn(elout(:,2)).^2).*(yn(elout(:,3))-...
    yn(elout(:,1)))+...
    (xn(elout(:,3)).^2+yn(elout(:,3)).^2).*(yn(elout(:,1))-...
    yn(elout(:,2)));
ycnum=(xn(elout(:,1)).^2+yn(elout(:,1)).^2).*(xn(elout(:,3))-...
    xn(elout(:,2)))+(xn(elout(:,2)).^2+...
    yn(elout(:,2)).^2).*(xn(elout(:,1))-xn(elout(:,3)))+...
    (xn(elout(:,3)).^2+yn(elout(:,3)).^2).*(xn(elout(:,2))...
    -xn(elout(:,1)));
elout(:,9)=xcnum./den/2;
elout(:,10)=ycnum./den/2;
elout(:,11)=(xn(elout(:,1))-elout(:,9)).^2+...
    (yn(elout(:,1))-elout(:,10)).^2;
elout(:,12)=mat;
s=2/3*((yn(elout(:,1))+yn(elout(:,3))-2*yn(elout(:,2)))...
    .*(xn(elout(:,2))-xn(elout(:,1)))-(xn(elout(:,1))+...
    xn(elout(:,3))-2*xn(elout(:,2))).*(yn(elout(:,2))-...
    yn(elout(:,1))))./den;
elout(:,13)=xn(elout(:,1)).*(1-s)+1/2*s.*(xn(elout(:,2))...
    +xn(elout(:,3)));
elout(:,14)=yn(elout(:,1)).*(1-s)+1/2*s.*(yn(elout(:,2))...
    +yn(elout(:,3)));
